﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace DynaPrecision.LinearSystem
{
	/// <summary>
	/// LU decomposition of a matrix
	/// </summary>
	public class DoolittleDecomposition
	{
		/// <summary>
		/// Decompose the decomposition, halfbandwidth is used
		/// </summary>
		/// <typeparam name="TInputMatrix"></typeparam>
		/// <typeparam name="TOutputMatrix"></typeparam>
		/// <typeparam name="TData"></typeparam>
		/// <param name="input"></param>
		/// <param name="result"></param>
		/// <param name="HalfBandwidth"></param>
		public static void Decompose<TData>(IMultiDimensionalArray<TData> input, IMultiDimensionalArray<TData> result, int HalfBandwidth)
			
		{
			IMultiDimensionalArray<TData> L = result, U = result;
			int i = 0, j = 0, k = 0;
			int totalI = input.GetLength(0), totolJ = input.GetLength(1);
			for (; i < totalI; i++)
			{
				j = i;
				//Update the Upper one
				for (; j < totolJ && j < i + HalfBandwidth + 1; j++)
				{
					U[i, j] = input[i, j];
					k = 0 > i - HalfBandwidth ? 0 : i - HalfBandwidth;
					for (; k < i; k++)
					{
						U[i, j] = U[i, j] - L[i, k] * (dynamic)U[k, j];
					}
				}
				//Update the Lower one
				j = i + 1;
				for (; j < totalI && j < i + HalfBandwidth + 1; j++)
				{
					L[j, i] = input[j, i];
					k = 0 > i - HalfBandwidth ? 0 : i - HalfBandwidth;
					for (; k < i; k++)
					{
						L[j, i] = L[j, i] - L[j, k] * (dynamic)U[k, i];
					}
					L[j, i] = L[j, i] / (dynamic)U[i, i];
				}
			}
		}
		/// <summary>
		/// General array decomposition
		/// </summary>
		/// <typeparam name="TInputMatrix"></typeparam>
		/// <typeparam name="TOutputMatrix"></typeparam>
		/// <typeparam name="TData"></typeparam>
		/// <param name="input"></param>
		/// <param name="result"></param>
		public static void Decompose<TData>(IMultiDimensionalArray<TData> input, IMultiDimensionalArray<TData> result)
		{
			IMultiDimensionalArray<TData> L = result, U = result;
			int i = 0, j = 0, k = 0;
			int totalI = input.GetLength(0), totolJ = input.GetLength(1);
			for (; i < totalI; i++)
			{
				j = i;
				//Update the Upper one
				for (; j < totolJ; j++)
				{
					U[i, j] = input[i, j];
					k = 0;
					for (; k < i; k++)
					{
						U[i, j] = U[i, j] - L[i, k] * (dynamic)U[k, j];
					}
				}
				//Update the Lower one
				j = i + 1;
				for (; j < totalI; j++)
				{
					L[j, i] = input[j, i];
					k = 0;
					for (; k < i; k++)
					{
						L[j, i] = L[j, i] - L[j, k] * (dynamic)U[k, i];
					}
					L[j, i] = L[j, i] / (dynamic)U[i, i];
				}
			}
		}

		public static void Decompose<TData>(IMultiDimensionalArray<TData> input, IMultiDimensionalArray<TData> result, int HalfBandwidth, int LowerFullWidthCount)
		{
			IMultiDimensionalArray<TData> L = result, U = result;
			int i = 0, j = 0, k = 0;
			int totalI = input.GetLength(0), totolJ = input.GetLength(1);
			for (; i < totalI - LowerFullWidthCount; i++)
			{
				j = i;
				//Update the Upper one
				for (; j < totolJ; j++)
				{
					if (!((j - i) < (HalfBandwidth + 1) || j + LowerFullWidthCount >= totolJ)) continue;
					U[i, j] = input[i, j];
					k = 0 > i - HalfBandwidth ? 0 : i - HalfBandwidth;
					for (; k < i; k++)
					{
						U[i, j] = U[i, j] - L[i, k] * (dynamic)U[k, j];
					}
				}
				//Update the Lower one
				j = i + 1;
				for (; j < totalI; j++)
				{
					if (!((j - i) < (HalfBandwidth + 1) || j + LowerFullWidthCount >= totolJ)) continue;
					L[j, i] = input[j, i];
					k = 0 > i - HalfBandwidth ? 0 : i - HalfBandwidth;
					for (; k < i; k++)
					{
						L[j, i] = L[j, i] - L[j, k] * (dynamic)U[k, i];
					}
					L[j, i] = L[j, i] / (dynamic)U[i, i];
				}
			}
			totalI = input.GetLength(0); totolJ = input.GetLength(1);
			for (; i < totalI; i++)
			{
				j = i;
				//Update the Upper one
				for (; j < totolJ; j++)
				{
					U[i, j] = input[i, j];
					k = 0;
					for (; k < i; k++)
					{
						U[i, j] = U[i, j] - L[i, k] * (dynamic)U[k, j];
					}
				}
				//Update the Lower one
				j = i + 1;
				for (; j < totalI; j++)
				{
					L[j, i] = input[j, i];
					k = 0;
					for (; k < i; k++)
					{
						L[j, i] = L[j, i] - L[j, k] * (dynamic)U[k, i];
					}
					L[j, i] = L[j, i] / (dynamic)U[i, i];
				}
			}
		}
	}
}